Hit the ground running: an introductory example ================================================ This chapter mirrors the canonical introductory example used in the legacy documentation. The same ``model.rs`` file and the same calibration / priors / sampling code are used; only the constructor (``dsge_model(...)`` instead of ``rise(...)``) and the named-argument form of ``set`` / ``estimate`` differ. Reading this chapter against the legacy *Teaser Example* makes the stable-to-modern syntax delta concrete in one place. .. contents:: :local: :depth: 2 The model economy ------------------ A three-equation New Keynesian model. The first equation is an Euler-type IS curve relating output to the real interest rate. .. math:: :label: IS curve \beta E_{t}\frac{r_{t}}{\pi _{t+1}}\left[ \frac{1-\frac{\kappa }{2}\left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}{1-\frac{\kappa }{2} \left( \pi _{t+1}-\pi _{t}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}\right] \frac{Y_{t}}{Y_{t+1}}=1 The second equation is a nonlinear Phillips curve with price adjustment costs. .. math:: :label: Phillips curve \begin{array}{ccc} 0 & = & \chi _{0}\eta _{t}\left[ 1-\frac{\kappa }{2}\left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) ^{2}\right] \left( \frac{Y_{t}}{Z_{t}}\right) ^{1+\chi } \\ & & +\left( 1-\eta _{t}\right) -\kappa \left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) \pi _{t} \\ & & +\beta \kappa E_{t}\left( \pi _{t+1}-\pi _{t}^{\alpha }\pi ^{1-\alpha }\right) \pi _{t+1}\left[ \frac{1-\frac{\kappa }{2}\left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}{1-\frac{\kappa }{2}\left( \pi _{t+1}-\pi _{t}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}\right] \end{array} The third equation is a monetary policy reaction function with a state-dependent inflation reaction. .. math:: :label: Taylor rule \frac{r_{t}}{r}=\left( \frac{r_{t-1}}{r}\right) ^{\rho }\left[ \left( \frac{ \pi _{t}}{\pi _{\ast }}\right) ^{\psi _{\pi }\left( s_{t}\right) }\right] ^{\left( 1-\rho \right) }\exp \left( \sigma _{r}\varepsilon _{r,t}\right) The fourth equation is the cost-push shock process. .. math:: :label: Cost-push shock process \log \eta _{t}=\left( 1-\rho _{\eta }\right) \log \eta +\rho _{\eta }\log \eta _{t-1}+\sigma _{\eta }\varepsilon _{\eta ,t} The fifth equation is the technology shock process. .. math:: :label: Technology shock process \log \left( \frac{Z_{t}}{Z_{t-1}}\right) =\left( 1-\rho _{z}\right) \mu _{z}\left( s_{t}\right) +\rho _{z}\log \left( \frac{Z_{t-1}}{Z_{t-2}}\right) +\sigma _{z}\varepsilon _{z,t} Regime switching enters through :math:`\psi_{\pi}` (the central bank's inflation reaction) and :math:`\mu_{z}` (the growth rate of technology), switching in lockstep on a chain we name ``snk``. Uppercase variables are non-stationary; lowercase variables are stationary. The model file --------------- The same ``model.rs`` file used in the legacy documentation. The file lives next to this chapter; full listing: .. literalinclude:: model.rs :language: rise The two distinguishing declarations: * ``@endogenous(log) Y "Output", R "interest rate", PAI "Inflation", Z "Technology", ETA "Cost-push shock process", DY "Output growth"`` -- every endogenous variable is log-linearised. * ``@parameters(snk,2) muz "growth-rate of technology", psi_pai "reaction function: inflation"`` -- the two parameters that switch, on a two-state chain named ``snk``. The full file is `model.rs `_. Building the model object -------------------------- :: m = dsge_model('model'); (The legacy documentation used ``m = rise('model')`` or ``m = dsge('model')``. The modern toolbox replaces the class-named constructor with shape-specific factories; ``dsge_model`` is the DSGE factory.) Without further options, RISE computes and stores first-order symbolic derivatives suitable for a first-order perturbation. Parameterizing the model ------------------------- :: p = struct(); p.beta = 0.99; p.kappa = 161; p.paitarg = 1.02^0.25; p.alpha = 0.5; p.eta = 6; p.chi = 0.7; p.chi0 = 1; p.psi_pai_snk_1 = 2.5; p.psi_pai_snk_2 = 0.9; p.rho = 0.7; p.rhoz = 0.75; p.rhoeta = 0.75; p.muz_snk_1 = 1.04^.25; p.muz_snk_2 = 1.01^.25; p.sigz = 0.05; p.sigeta = 0.05; p.sigr = 0.05; p.snk_tp_1_2 = 1 - 0.95; p.snk_tp_2_1 = 1 - 0.9; m = set(m, parameters = p); Switching parameters carry the chain name and state in their name: ``psi_pai_snk_1`` is :math:`\psi_{\pi}` in state 1 of chain ``snk``; ``snk_tp_1_2`` is the transition probability from state 1 to state 2 of chain ``snk``. Solving the non-stationary model --------------------------------- The model is non-stationary; we let RISE solve it on the balanced-growth path with the Maih-Waggoner perturbation strategy and the matching Newton solver:: m = solve(m, ... solve_perturbation_type = 'mw', ... solve_bgp = true, ... solver = 'mn'); print_solution(m); A one-shot ``stoch_simul`` run ------------------------------- :: info = stoch_simul(m); ``info`` carries simulated data, moments, IRFs, autocorrelations, skewness and kurtosis. Impulse responses ------------------ :: myirfs = irf(m); quick_irfs(m, myirfs, get(m, 'endo_list(original)')); The output is a struct of ``ts`` keyed first by shock name then by endogenous variable; ``endo_list(original)`` excludes the auxiliary variables RISE introduces internally. 5th-order perturbation ----------------------- :: m5 = solve(m, ... solve_order = 5, ... solve_derivatives_type = 'automatic'); print_solution(m5); The model was originally parsed for first-order symbolic derivatives; for fifth order we switch the derivative engine to automatic (algorithmic) differentiation. Alternatively, re-parse the model with ``max_deriv_order = 5`` to get symbolic derivatives up to order 5. Data from FRED --------------- :: xrange = '1960Q1:2022Q3'; rawdb = fetch_fred({'BPCCRO1Q156NBEA', ... 'BOGZ1FL072052006Q', ... 'NAEXKP01USQ657S'}); db = struct(); db.PAI = 1 + rawdb(1).series(xrange)/100; db.R = 1 + rawdb(2).series(xrange)/100; db.GROWTH = 1 + rawdb(3).series(xrange)/100; These are the canonical PCE-inflation, federal funds rate, and US-output series used in the documentation example. Priors ------- The legacy chapter sets priors in **quantile** form. We do the same. Several priors illustrate the four-, five-, and six-slot variants (see :doc:`../Estimation/Main Estimation` for the full parametrisation taxonomy):: priors = struct(); % 4-slot quantile form priors.kappa = {p.kappa, 5, 20, 'gamma(.9)'}; priors.alpha = {p.alpha, 0.05, 0.948,'beta(.9)'}; priors.chi = {p.chi, 1.5, 3, 'gamma(.9)'}; % 6-slot generalized-beta form (extra hyperparameters) priors.eta = {p.eta, 3, 8, 1, 12, 'beta(.9)'}; priors.muz_snk_1 = {p.muz_snk_1, 1.02^.25, 1.05^.25, 1, 1.07^.25, 'beta(.9)'}; priors.muz_snk_2 = {p.muz_snk_2, 1.01^.25, 1.03^.25, 1, 1.04^.25, 'beta(.9)'}; % 4-slot quantile form priors.rho = {p.rho, 0.05, 0.948, 'beta(.9)'}; priors.rhoz = {p.rhoz, 0.05, 0.948, 'beta(.9)'}; priors.rhoeta = {p.rhoeta, 0.05, 0.948, 'beta(.9)'}; % 4-slot + hard truncation (6-slot) priors.sigz = {p.sigz, 0.0005, 1.0, 'sichisq(.9)', 0, 3}; priors.sigeta = {p.sigeta, 0.0005, 1.0, 'sichisq(.9)', 0, 3}; priors.sigr = {p.sigr, 0.0005, 1.0, 'sichisq(.9)', 0, 3}; % 4-slot, transition probabilities priors.snk_tp_1_2 = {p.snk_tp_1_2, 0.01, 0.411, 'beta(.9)'}; priors.snk_tp_2_1 = {p.snk_tp_2_1, 0.01, 0.411, 'beta(.9)'}; % 6-slot generalized-beta, switching reaction priors.psi_pai_snk_1 = {p.psi_pai_snk_1, 1.1, 2.5, 1, 3, 'beta(.9)'}; priors.psi_pai_snk_2 = {p.psi_pai_snk_2, 0.5, 1.1, 0.1, 1.5, 'beta(.9)'}; Visualizing the priors ----------------------- :: plotOpts = struct(); plotOpts.prior_trunc = 2.6e-3; rdist.plot(priors, plotOpts, [], struct('linewidth', 2)); Posterior maximisation ----------------------- :: mest = estimate(m, ... data = db, ... estim_priors = priors, ... kf_init_variance = 1, ... kf_presample = 10, ... optimizer = 'bee_gate'); * ``kf_init_variance = 1`` -- initial variance of every endogenous variable, set to 1 because the model is non-stationary and does not have a finite unconditional variance. * ``kf_presample = 10`` -- the first 10 observations are not included in the likelihood; this is a defensive complement to ``kf_init_variance``. * ``optimizer = 'bee_gate'`` -- the artificial bee colony global-search optimizer, robust on this calibration. Posterior simulation --------------------- We use the random-walk Metropolis-Hastings sampler from the modern ``+rsamplers`` package. ``pull_objective`` returns a value to **minimize**; ``rsamplers.rwmh`` takes a value to **maximize**, so we negate the objective before passing it in:: [objective, lb, ub, x0, SIG] = pull_objective(mest); scale = 0.0794; myOpts = struct(); myOpts.tunedCov = scale * SIG; myOpts.N = 20000; energy = @(varargin) -objective(varargin{:}); results = sample(rsamplers.rwmh(energy, x0, lb, ub, myOpts)); results{1}.stats Marginal data density via the bridge estimator ----------------------------------------------- RISE ships nine ways of computing the marginal data density. For the Meng-Wong (1996) bridge estimator:: mddobj = mdd(results, energy, lb, ub, [], [], true); bridge(mddobj, true, mdd.global_options); The full driver ---------------- The complete driver script is `model_driver.m `_ shipped alongside this chapter. Running it end to end requires network access for the FRED fetch. The legacy documentation ships the same model file and a driver with the same content under the legacy syntax. Diffing the two drivers side by side is the cleanest way to see the stable-to-modern syntax delta in practice.